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Abstract 


A low-energy neutron transport, algorithm for use in space radi- 
ation protection is developed. The algorithm is based upon a multi- 
group analysis of the straight-ahead Boltzmann equation by using a 
mean value theorem for integrals. This analysis is accomplished by 
solving a realistic but simplified neutron transport test problem. The 
test problem is analyzed by using numerical and analytical proce- 
dures to obtain an accurate solution within specified error bounds. 
Results from the test problem arc then used for determining mean 
values associated with rescattering terms that are associated with a 
multigroup solution of the straight-ahead Boltzmann equation. The 
algorithm is then coupled to the Langley HZETRN code through the 
evaporation source term. Evaluation of the neutron flue nee gener- 
ated by the solar particle event of February 53 } 1956, for a water and 
an aluminum- water shield-target configuration is then compared with 
LAHET and MCNPX Monte Carlo code calculations for the same 
shield-target configuration. The algorithm developed showed a great 
improvement in results over the unmodified HZETRN solution. In 
addition , a two- directional solution of the evaporation source showed 
even further improvement of the fluence near the front of the water 
t a rge t wh e re d iff us i o n fro m t h e fw nt surface is impo rt a n t . 


Introduction 

The purpose of this paper is to present an improved algorithm for the analysis of the transport 
of low-energy neutrons arising in space radiation protection studies. The design and operational 
processes in space radiation shielding and protection require higlily efficient computational 
procedures to adequately characterize time-dependent environments, time-dependent geometric 
factors, and to address shield evaluation issues in a multidisciplinary' integrated engineering 
design environment. One example is the recent study of the biological response in exposures 
to space solar particle events (SPE’s) in which the changing quality of the radiation fields at 
specific tissue sites is followed over 50 hours of satellite data to evaluate time-dependent factors 
in biological response of the hematopoietic system (ref. 1). Similarly, the study of cellular 
repair dependent effects on the neoplastic cell transformation of a C3H10T ^ population in low 
Earth orbit, where trapped radiations and galactic cosmic rays vary continuously in intensity and 
spectral content about the orbital path (ref. 2), requires computationally efficient codes to match 
time-dependent boundary conditions around the orbital path. But even in a steady environment 
which is homogeneous and isotropic, the radiation fields within a spacecraft have large spatial 
gradients and highly anisotropic factors so that the mapping of the radiation fields within the 
astronaut's tissues depends on the astronaut timeline of location and orientation within the 
spacecraft ulterior where large differences in exposure patterns that depend on the activity of 
the astronaut have been found (ref. 3). Obvious cases exist where rapid evaluation of exposure 
fields of specific tissues are required to describe the effects of variations in the time-dependent 
exterior environment or changing geometric arrangement. A recent study of the time-dependent 
response factors for 50 hours of exposure to the SPE of August 4, 1972, required 18 CPU hours 
on a VAX 4000/500 computer by using the nucleon-light ion section of the deterministic high 
charge and energy transport code HZETRN. The related calculation with a standard Monte 
Carlo code such as HETC or LAHET, which only handles neutrons, protons, pions, and alphas, 
would have required approximately 2 years of computer time to complete the study. The design 
environment also requires rapid evaluation of the radiation fields to adequately determine effects 
of multiparameter design changes on system performance (refs. 4 and 5). These effects are the 


driving factors in the development and use of deterministic codes and in particular the HZETRN 
code system that handles all naturally occurring atomic ions and neutrons. 

The basic philosophy for the development of the deterministic HZETRN code began with the 
study by Alsmiller et al. (ref. 6) with an early version of HETC, wherein they demonstrated 
that the straight-ahead approximation for broad beam exposures was adequate for evaluation of 
exposure quantities. Wilson and Khandelwal (ref. 7) examined the effects of beam divergence 
on the estimation of exposure in arbitrary convex geometries and demonstrated, that the errors 
in the straight- ahead approximation are proportional t-o the square of the ratio of the beam 
divergence to the radius of curvature, which is small in typical space applications. From a 
shielding perspective, the straight-ahead approximation overestimates the transmitted flux, and 
the error is found to be small in space radiation exposure quantities. Our first, implementation 
of a numerical procedure w T a,s performed by Wilson and Lamkin (ref. 8) as a numerical iterative 
procedure of the charged components perturbation series expansion of the Boltzmann transport 
equation and showed good agreement with Monte Carlo calculations for modest penetrations to 
where neutrons play an important role. The neutron component was added by Lamkin (ref. 9); 
this closed the gap between the deterministic code and the Monte Carlo code. The resulting 
code was fast compared with the Monte Carlo codes but still lacked efficiency in generating and 
handling large data arrays, which would be solved in the next generation of codes. 

The transport of high-energy ions is well adapted to the straight-ahead approximation. In 
fact, a more common assumption that secondary ion fragments are produced with the same 
velocity as the primary initial ion (ref. 10) is inferior to the straight-ahead approximation 
contrary to intuition (ref. 11). The Boltzmann transport equation for the particle fields <j>j(x^ E) 
is given in the straight-ahead and continuous slowing down approximations as 


d_ 

dx 


Sj(E) + cTj(E)]<p j (x,E)= r q jk (E,E?)<t> k { x,E') dE' 


( 1 ) 


where x is the depth of penetration, E is the particle kinetic energy, Sj(E) is the particle 
stopping power, &j{E) is the macroscopic interaction cross section, and aj^E^E*) is the 
macroscopic cross section for particle k of energy E f produced as a result of the interaction 
with a particle j of energy E . At Langley Research Center for all the code development, it 
has been customary to invert the differential operator and implement it exactly as a marching 
procedure (ref. 12), and the remaining issue has been in approximating the integral term on the 
right-hand side of equation (1). The implementation for the heavy fragments was facilitated by 
the assumption that the fragment velocity is the same as the primary ion which is inadequate 
for the description of the coupled nucleonic and light ion components. A compatible nucleonic 
transport procedure was developed by Wilson et al. (ref. 13) and showed good agreement with 
exposure quantities evaluated by Monte Carlo transport procedures (ref. 14). The transport of 
the nucleonic component was developed by assuming that the midpoint energy within the step 
was the appropriate energy to evaluate the integral term. Thus, the residual range of the proton 
will reduce by hj 2 before the interaction and the secondary proton residual range will reduce 
by h/2 before arriving at the next marching step. Neutrons show no loss in residual range as 
their stopping power is zero. This choice was shown to minimize the second-order corrections to 
the marching procedure (ref. 15). Although reasonable agreement on exposure quantities from 
Monte Carlo calculations was obtained, the resultant neutron flux at the lowest energies was 
substantially below the Monte Carlo result in the range of 0.01 to several MeV and required 
improvement (ref. 16). Analysis concluded that the problem was in the rescattering terms in 
wliich the number of elastic scattered neutrons was underestimated numerically, which must be 
addressed as suggested by Shinn et al. (ref. 16). 
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The issue of evaluation of the integral term of the Boltzmann equation for the elastic scattering 
is the next issue to be resolved in the development of the HZETRN deterministic code. Once 
the elastic scattering events are adequately represented and the associated improvements in the 
neutron flux are made, one still needs to address the issue of the adequacy of the nuclear database 
for nucleonic transport in the HZETRN code system (ref. 13). 

Formulation of Transport Equations 

Define the differential operator B as 

= + „ j(E) 4(x,E) (2) 

and consider the following one- dimensional Boltzmann equation from reference 17 

b \h ^ n r °jk( E £ } ) M*#’) dE * (3) 

k 

where <f> j is the differential flux spectrum for the type j particles, Sj{E) is the stopping power 
of the type j particles, and <rj(E) is the total macroscopic cross section. The term <rjk(E,E?), 
a macroscopic differential energy cross section for redistribution of particle type and energy, is 
written as 

a jk (E,E') = Pi) fjkd E ’ E ') 

13 

where fjk^{E,E f ) is the spectral redistribution, aj is a microscopic cross section, and ps is the 
number density of {3 type atoms per unit mass. The spectral terms are expressed as 

fjk,t. ) = fj\,3 + fjkjl + fjk ti 3 

where /jj. j represents the elastic redistribution in energy, /j'^ ^ represents evaporation terms, 
and fjj. ,j represents direct knockout terms. The elastic term is generally limited to a small 
energy range near that of the primary particle. The evaporation process dominates over the low 
energies (E < 25MeV) and the direct cascading effect dominates over the high energy range 
(E > 25 MeV) as illustrated in figure 1. 

Equation (3) is then written for j = n as 

BM = E r E (V + r.k,a + lik,a) *»(*.«') i* w 

t JE ll 

wliich is expanded to the fonu 


» m = r e p* ,,e ' 

JE p 

+ XI / + ft,k,J + fnk,8> < l > k( J ’^’) dE' 

k£n JE (i 


(5) 
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Define the integral operators I as 


’w = r E w *?&) nlj dE ' 

JE p 

li k) [*] = E P0 fnk,3 <K*&) dE ' 

J E p 

i l d k) w = r E pn *?&) fit# tW) dE ' 

J E a 


where k = n denotes coupling to neutron collisions and k — p denotes the neutron source from 
proton collisions. When considering only neutrons and protons, equation (5) can be written in 
the linear operator form as 


b[« = 4"'w + /n*»] + /;■' w + /;:>,] + r'w + rw 


r(«) 


(p)r 


At>) r 


r(p)r 


( 6 ) 


Note that /^[<^,] does not contribute to the neutron field; therefore, equation (6), with <j> n 
replaced by <£, is written as 

BM = ^"V] + 4" ) M + /J" ) W + rf'' ) W + 4 ,,| W ( 7 ) 

Assume a solution to equation (7) of the form <f> — <f> t where <p t is the solution for evaporation 
sources and contributes over the low-energy range and is the solution for the direct knockout 
sources and contributes mainly over the high-energy range as suggested by figure 1. Substitute 
this assumed solution into equation (7) and find 

Bw = sw + sw = 4;‘ Vi + ’w + d"Vi + il"\u 

+ ;<">[#,] + + /<">[*,] (8) 


The terms l\ U \<p t ] and I^ l \<Pt] are near zero and are ignored because evaporation neutrons at 
low energies do not produce additional evaporation neutrons, and the direct cascade effects have 
very small cross sections over the low-energy range of <f> t and hence does not contribute any 
production over the low- or high-energy range, Further assume that <p ( j is calculated by the 
TTZETRN program so that <f) t { is a solution of the equation 

B[<t>d] = + I { /%] ( 9 ) 

Tliis assumption simplifies equation (8) to the form 

B[* e ] = £\*'] + il n) m + A ,,] y> P ] (io) 

Define the elastic scattering terms 

= Pfi <r,j(E') fjl,j,{E,E') 


with units of cm 2 /g-MeV, and note that for neutrons the stopping power Sj (E) is zero and 
equation (10) reduces to the in tegro- differential transport equation with source term 


‘ d 

dx 


+ v(E) 



<t s , 3 (E,E') M*X') dE ’ + g(E,x) 


( 11 ) 
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Equation (11) represents the steady-state low-energy neutron fluence <f) t (x,E) at depth x and 
energy E . The various terms in equation (11) are energy E with units of MeV, depth in medium 
is x with units of g/cm 2 , <f> t (x } E) (in particles/cm 2 -MeV) is the evaporation neutron fluence, 

and g(E, x) - + (in particles/g-MeV) is a volume source term to be evaluated by 

the HZETRN algorithm. Equation (11) is further reduced by considering the neutron energies 
before and after a collision. The neutron energy E n after an elastic collision with a nucleus of 
mass number Ap^ initially at rest, is, from reference 18, 


En = E 


Aj* T cos 9 + 1 

( A Tp + *) 2 


( 12 ) 


where E is the neutron energy before the collision, A^. is the atomic weight of the ith type of 
atom being bombarded, and 9 is the angle of scatt er. Define the ratio 


Q d = 



(13) 


as a constant less than 1 and note that when 0 = 0, E n = E y and when 0 — tt, E n = Ea j. 
Therefore, change the limits of integration in equation (11) to ( E , E/aj) which represent the 
kinetically allowed energies for the scattered neutron to result in an energy E. Equation (11) 
then is written as 


d__ 

dx 


+ <?(E) 


4>i 


(^) = E/ £ 




<r^(E,E!) 4> e (x,Ef) dE’ + g(E,x) 


(14) 


The quantity o in cnr / g is a macroscopic cross section given by 

<r=Ew ( 15 ) 

d 

where p,j is the number of atoms per gram and a ^ is a microscopic cross sect ion in cnr/atom. 
Reference 19 provides approximate Maxwellian averages of cross-section values in barns which 
are used herein for studies of solution techniques. These values are listed in table 1 along with 
other parameters of interest for selected elements. Other units for equation (11) are obtained 
from the previous units by using the scale factor representing the density of the material in units 
of g/crn 3 . 


Table 1. Parameter Values for Selected Elements 


Element 

A T ;) 

Elastic 

cross section, 
barns" 

Density, 
g/ cm 3 

a H 

Lithium, Li 

7 

1.050 

0.534 

0.563 

Carbon, C 

12 

4.739 

0.352 

0.716 

Aluminum, A1 

27 

1.348 

2.7 

0.862 

Calcium, Ca 

40 

2.99 

1.54 

0.905 

Iron, Fe 

56 

11.40 

7.85 

0.931 

Lead, Pb 

207 

11.194 

11.342 

0.981 


"Maxwellian averages (ref. 19). 


Mean Value Theorem 


Throughout the remaining discussions, the following mean value theorem is used for integrals. 

Mean Value Theorem: For <p(x,E) and f{E) continuous over an interval a < E < b such 
that (1) <p(xji) does not change sign over the interval (a, 6), (2) <p[x y E) is integrable over the 
interval (a, b), and (3) f(E) is bounded over the interval (a, 6), there exists at least one point c 
such that ^ ^ 

' [■ f(E) <j>(x,E) dE = /(e) f <p(x y E) dE (a<e< b) 

Ja J a 

In particle transport, this mean value approach is not commonly used. In reactor neutron 
calculations, an assumed spectral dependence for ^x^E) is used to approximate the integral over 
energy groups. The present use of the mean value theorem is free of these assumptions; thus, 
more flexibility is allowed in the HZETRN code, and the result is a fast and efficient algorithm 
for low neutron analysis. 


Multigroup Method 

Consider the case where there is only one value of (3 wliich represents neutron penetration 
into a single element material and let <f>t be denoted by <j>. Equation (14) is integrated from E % 
to Ei+ j with respect to the energy E to obtain 



+ !«*(*,£) 

dx 


f&i + 1 

dE + J (t{E) <f>(x,E) 


dE=I i +i i 


( 16 ) 


where 


I; = 



<r,jj{E t Ef) <p(x,E') dE' dE 


(17) 


and 



(18) 


As a test case for developing solution techniques, we use the approximate source and scattering 
terms taken from subroutine FBEKT of the HZETRN code (ref. 5), g = g(E,x) = KEe 
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with K and T constants, and the elastic scattering term from subroutine ELSPEC of the 
HZETRN code (ref. 5), 

ViAW) = “ 

with r constant, so that equation (18) Is easily integrated to obtain 


(E 1, )re~ r ^'~ g) 
1 _ e (l — «)r£' 


= KT (Eie~ Ei > T - E i+ \ + KT 2 (c _E ’/ T - e _E '+ l/r ) 


(19) 


The quantity 

$;(*) = f ’ ,+1 Hx,E) dE (20) 

J Ej 

is associated with the i th energy group, so that — — — d>,(a:) represents an average flue nee for 

fl E'i 

each energy group. Then equation (16) can be written in temis of T, (.r) as follows. In the first 
term of equation (16), interchange the order of integration and differentiation to obtain 


f E >+ 1 d<j)(x,E) 

L 


dE = 


dd>j(x) 

dx 


( 21 ) 


With the previously stated mean value theorem for. integrals, the second term in equation (16) 
can be expressed as 

f E i+i 

/ <r4>(x£)dE=a*;(x) (22) 

JEi 

where JF = <r[Ej. + 0(Ej + \ — E;)], for some value of 6 between 0 and 1. 

For the term I; in equation (17), the order of integration Is interchanged. Various partitioning 
schemes are illustrated in figure 2. The integration of equation (17) depends upon the energy 
partition selected. For example, figure 2(b) illustrates an energy partition where Ej+\ < EJa , 
and in this case, equation (17) can be written as 


/,•= / £,+l f* 

JE'=E t JE=Ej 


fEj/a pE i+l ,E i+ [/« pE i+l 

HdEdE’+ / / HdEdE'+ / II dE dE' (23) 

JE'=Ei + \ JE=Ei JE’=E,/o Je=uE’ 


where H = a s (E,E') <f>(x,E'). Figure 2(c) depicts the case where E i+ j = EJa exactly for all i. 
In this special case, equation (17) reduces to 


Ii 


f E i + 1 f E . f E i+ 1/« f E <+ 1 , 

/ / H dE dE + / / H dEdEf 

Je'=E, JE=Ei JE'=E, + \ J E=uE' 


(24) 


The selection of an energy partition can lead to two or more distinct groups associated with 
each interchange in the order of integration (for example, see fig. 3). The integrand II can be 
integrated with respect to E and the results expressed in terms of the quantities 



— e 


TO 


and 


G(E') = 


ojE) e~ rE ' 

1 _ e -(l-a)rE' 
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and equation (24) can be written in the form 


7; - f E,+] G(E') F{E\Ei ) <p{x,E') dE' 

J E'=Ej 

+ f E,+]/ " G(E’) F(E i+ uaE') <j>(x,E r ) dE' (25) 

Je'=e, +] 

To illustrate the basic idea behind the multigroup method, use the same mean value theorem 
for integrals and write equation (25) in the form 

Ii=G(E*) F(E* ,Ej) + G(E* +l ) F(E l+l ,aE* +l ) 4> ;+1 

where Ej < E* < E;/a and < E*+± < E;_ The special partitioning of the energy as 
illustrated in figure 2(c) enables us to obtain from equation (16) a system of ordinary differential 
equations as follows: 



$0 


"«11 «12 


' <*>0 " 


& 

d 

$1 


o-22 023 -0- 








— 




+ 


dx 

$.¥-2 


-0- flJV— 1,.V— 1 O.V-1..V 


< t > jY-2 

£.V-2 



-$JV— 1- 


a N N - 




-t.w-l - 


(26) 


where ajj — G(E *) F(E* ,Ej) — 7r and Qij+i ~ G ( E* + j ) F ( Ej+ i , <*Ej+ j ) . Further assume that 
for large values of N , 4>; = 0 for all i > N. This assumption gives rise to the following system of 
ordinary differential equations: 


dV 

dx 


Ay + b 


subject to the initial conditions y( 0) = Here y is the column vector of 4>; values, 
col($ 0 ^l»-- ■ s ^ A r — 1 ), the matrix A is an N by N upper triangular matrix, and b is the column 
vector col (£q,4i, . . . l). Hi a similar manner, the integrals in equation (23) can be evaluated 

for other kinds of energy partitioning and a system of equations having tin" form of equation 
(26) obtained. However, for these other energy partitions, the structure of the N by N square 
matrix A will change. It remains upper triangular but with more off-diagonal elements which 
depend upon the type of energy partition. (See, for example, fig. 3.) For our purposes the system 
of equations (eq. (26)) is used to discuss some of the problems associated with the multigroup 
method. 


Of prime concern is how an energy grid is to be constructed and how this energy grid controls 
the size of the matrix in equation (26). Consider the construction of the energy partition 


{ 


Ea 



J J 

a 


E o 

a 1 ' 


E o_ 

’a* 


} 


where Eo = 0.1 MeV, for the selected elements of lithium, aluminum, and lead. Table 2 
illustrates integer values of N necessary to achieve energies greater than 30 MeV. These values 
of N represent the size of the matrix associated with the number of energy groups. The value 
Eo = 0.1 MeV, in terms of human exposure, represents a lower bound where lower energies are 
not important . The value of 30 MeV represents an upper limit for the evaporation particles. 
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Table 2. Energy Partition Size N 


Ele ment 

Q 

N 

0.1 fa N 

Lithium 

0.563 

10 

31.53 

Aluminum 

0.862 

39 

32.75 

Lead 

0.981 

298 

30.38 


Observe that for energy partitions where E}+\ < E;/a the values of N are larger, and if 
£ 7+1 > Ei/ a the values of N arc smaller. The cases where Ej+ \ > Ej/o give rise to problems 
associated with the integration over the areas A\ and of figure 2 (d) when the order of 
integration is interchanged. In this figure, the area A\ is associated with the integral defining 
$•, and the area A 2 is a remaining area associated with an integral that is some fraction of the 
integral defining $;-fi which is outside the range of integration. Therefore, some approximation 
mast be made to define this fractional part. Tins type of partitioning produces errors, due to any 
approximations, but it has the advantage of greatly reducing the size of the N by N matrix A 
at the cost of introducing errors into the system of equations. A more detailed analysis of the 
energy partition can be found in reference 20 . 

The case of neutron penetration into a composite material gives rise to the case where j3 > 1 
in equation (11). In this special case, equation (17) becomes 

h = £ J^' +] P <r*j{E,Ef) <j>(x,E') dE' d.E 


Select a = max(ai,or 2 i *• and construct the energy partition where E;+ j = E\f a. Then 

obtain a system of differential equations having the upper triangular form: 



$0 


"oil 

a 12 

a 13 ' 

• a IN' 


$0 


' & ' 


*1 



a 22 

«23 * 

' a 2N 


*1 



d 







+ 






-0- 

«33 ’ 

a NN- 


_$V— 1 . 

7 

Uy 

1 




(27) 


Observe that for some arbitrary energy grouping we have, for the element hydrogen, a case 
where the value of oj is zero. In this situation we must integrate over many energy groups as 
illustrated in figure 3. Some type of approximations must be made when the order of integration 
Is interchanged, depending upon the selected energy partitioning. Also the problem of selecting 
the mean values associated with each of these integrations exists. 


Mean Value Determination 

Consider the case of neutron fluence in a single shield material with the energy partitioning as 
illustrated in figure 2(c). This case is where successive energy values are given by Z?; + i — Eja 
for all values of the index i as it ranges from 0 to N . Select a finite value for N large enough 
that the assumption <J>,v 0 holds true. The system of equations in equation (26) is then 

a closed system and we can solve for the last term and then march backwards to solve 
$ y_p 3 \v_2i . • . . 
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The nonzero elements ajj for matrix A in equation (26) consists of the diagonal elements 
and the first diagonal above the main diagonal. This gives the values 


fl;; =G(E*) F(E*,Ej) — W 

l =G(£J +1 ) F(Ei+ l} aE? +l ) 


for i = 1 } . . A T , where E* and E* +1 are selected mean values associated with the lower and 
upper triangles illustrated in figure 2(c). These mean values vary with energy and were selected 
so that the multigroup solution agrees with the numerical solution of the test problem. The 
values determined empirically were 


where 


and 


where 


Ef^Ei + BiiEi+i-E;) 

£*+1 = E i+] + 02(£|+2 “ ^i+l) 


f li + m n (E- En) ~ ^1 

(£■ > £7n) 

h = S 71 + 

m li(E — En) - <*1 

(£22 < E < Eu) 

l 73 + 

?7Ji3(£’- E22) ~ h 

(E < E22) 

( 72+ ™2l(E — Eh) 

(E>E 11) 

$2 = S 12 + rn22(E - E\ \ ) 

(E22 < E < Eu) 

l 74 + E-n) 

(E < £22) 

74 = 0.93 

m n = 0.0030485 

777 21 = 0.004355 

7-2 = 0.90 

mis = 0.2490258 

777*22 — 0.249026 

73 = 0.30 

77743 = -0.3937186 

777*23 = —0.255920 

74 = 0.27 

E u = 3.037829 

E 22 = 0.5079704 


and 6\ is 0.0 for lead, 0.02 for aluminum, and 0.075 for lithium. These values of 0 for the mean 
value theorems were determined by trial and error so that the multigroup curves would have the 
correct shape and agree with the numerical solution. These selections for the mean values are 
not unique. 


Solution Method in Shield Materials 

Consider the energy partition £7+1 = E[j a and the resulting system of differential equat ions 
(eq. (26)). The solution of this system of equations is obtained by first solving the last equation 
of the system. This equation has the form 

= a NN * N _ 1 +£.y— 1(*) (*v_i(0) = 0) 


and has the solution 


*V-l(*) = kv-i(0) + r^(s) e-*™' ds 

. J 0 J 

which implies 

rrQ+Ad‘ 

$.V-l(*0 + A*) = e aNN $y_i (*o) + e a ' v ' v( * 0+Al ' ) / £v-l(s) c“ a ' v ' v * ds 

J xn 
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Now consider each of the remaining equations above the last equat ion in equation (26). A typical 
equation from this stack has the form 


d$j - 1 
dx 


- aa<bi - 1 + fi(x) ($i- 1(0) = 0) 


(28) 


where fi(x) — £i( x )+ a i,i+l*bi{ x ) known, since ^/(tf) is calculated before <F/_i(.r). This typical 
equation has the solution 


- € aii J 


*;-i(0) + r f is) 
J 0 


which implies 

$;_l(x 0 + Ajp) = e a u** $j_i(* 0 ) + e a "' 


^ /„ \ , ,o,;(J'u+Aj') j J 0 + ^ c ~ a ii s 


r 


Observe that for the system of equations in equation (27), the solution technique is essentially 

the same with the exception that the right-hand side of equation (28) is replaced by a summation 

N 

of the previously calculated terms, so that fi(x) = £/(#) + g 'j 

j=i+ J 


Numerical Solution 

The solutions obtained from the system of equations (eq. (26) or (27)) depend upon the 
selection of mean values associated with each energy interval. The selection of these mean 
values is determined by examining the numerical solution in certain special cases. We obtain a 
numerical solution of equation (1 1) in the special case given by 


g = g(E,x) = KEe~ E l T 

where K (particles/cm 3 -MeV) and T (MeV) are constants. We construct the solution over the 
spatial domain x > 0 and energy range 0.1 < E < 80 MeV. This domain is discretized by 
constructing a set of grid points Xf ~ i Ax and Ej — j AE for some grid spacing defined by 
Ax and AE values being used. For ij integers, define ujj = <j)(xj< Ej), then the transport 
differential-integral equation (11) can be written in a discrete form as follows, with the starting 
values wq j = 0 and j = 0 being used. For the first step in Aj*, approximate the flux by the 
accumulation of the source over the first interval as 

u\j = AxKEj e~ E J /T (29) 

followed by the numerical calculation of the rescattering term 


w: 


E j I" <t{E’)t e ~ T(E '- E J ) 


1 _ c -(l -«)tE> 


u(xi,E r ) dE 1 


(30) 


for 7—1. After this first, and each successive step, integrals of the type v;j given by equation (30) 
are evaluated with Simpson's one-third rule. Evaluate equation (30) for all energies j = 0, 1 , . . 
and then use a two-step algorithm in a repetitive fashion to advance the solution. For values of a 
near 1, the numerical solution of equation (11) requires that AE become small. The low-energy 
spectrum then becomes difficult to calculate without special procedures, as cited in reference 17. 
In this case, a two-step modified Euler predictor-corrector scheme is used (refs. 21 and 22), which 
Is defined by 
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Second step: 


— E 

fl,j= V\ j+Eje j - (tu i j 

. _ / « 1 J + Ax /i j 0’ = 0) 

Ti2, ' y l 5 ( U 1J-1 + 7i lj+l) + Ax/jj ( j > 0) 

Tliird step: 

hj = V 2J + Ej e~ E J - <TU 2 j } 

«3 J= wij + 2Ax/ 2 j J 


(31a) 


(31b) 


Tlif second step is an adoption of the Fredrichs method from reference 21. The third step is a 
central difference second-order st ep in Ax. After 100 applications of this two-step algorithm, we 
apply the following stability correction as suggested in reference 22: 


hj - l %j + E j e J ~ au zj 
«3J = \ ( W3 J + “2 j) + Ax hj 


(32) 


Note equations (32) are to be understood in an iterative sense and not strictly algebraic sense. 


Recursive Solution 

In the special case g(E,x ) = g(E), a solution to equation equation (11) is assumed of the 
form 

oo 

4>{x,E) = ME) /»(*) = 4>\ (E) f t (x) + <j> 2 (E) f 2 (x) +■■■ (33) 

n= 1 

Substitute this series into equation (11) and obtain a solution by requiring that <f> and / satisfy 


ME) = 9(E) 

K+l (E) = r ! ° h(E,E') ME') dE' 

f[(x) + afl(x) = 1 
f'r,(x) + cr/„(x) = /„_l(x) 


(34) 


for n = 1,2,3,..., where the differential equations are subject to the initial condition that 
/„(0) = 0 for all n. Here the terms for <f>n(E) are defined recursively and take a great deal of 
computational time for large values of n. The differential equations have the solutions given by 
the recursive relations 


/,(*)=!( l-e“") 

/„(>■)= r h- 1<“) 

J 0 


(35) 


which are easily evaluated for as large a value of n as desired. We find numerically that \fn{%)\ 
decreases with increasing n for x < 1 and increases for x > 1 so that the series solution does 
not converge in tills case. For \x\ < 1, we calculated the solution given by equation (35) for 
terms through n - 5 and n = 6 and compared them with the numerical solution. The mean 
values associated with the multigroup method were then adjusted so that the multigroup method 
agreed with the numerical solution and recursive solution for this special test problem. We then 
used these same mean values which where associated with numerical source terms as provided 
by the HZETRN code. 
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Comparison of Multigroup and Other Solutions 

The numerical solutions and recursive solutions of the test problem were then compared with 
the multigroup solution for neutron penetration in lithium, aluminum, and lead mediums. The 
results are illustrated in figures 4, 5, and 6. Excellent agreement is obtained in these three 
cases. In these figures, the solid line represents the numerical solution. The circles represent the 
recursive solution and the triangles represent the multigroup solution. The various curves were 
calculated for depths x of 0.1, 0.5, 1 .0, 5.0, 1CK0, 50.0 and 100.0 g/crn 2 . 

The multigroup method has huge advantage in its very short computational time needed to 
calculate the solution without loss of accuracy. The multigroup method takes less than 1 min 
of computational time, whereas the Monte Carlo methods require many hours of computa- 
tional time. 


Application for AI-H2O Shield-Target Configuration 

The previous development is now applied to an application of the multi group method 
associated with an aluminum- water shield-target configuration. In particular, consider the case 
where the source term g{E y x) in equation (11) represents evaporation neutrons produced per unit 
mass per MeV and is specified as a numerical array of values corresponding to various shield- 
target thicknesses and energies. The numerical array of values is produced by the radiation 
code HZETRN developed by Wilson et al. (ref. 23). The numerical array of values are actually 
given in the form g{E l ,x J ,y^] in units of particles/g-MeV, where y^ represents discrete values 

for various target thicknesses of water in g/cm 2 , xj represents discrete values for various shield 
thicknesses of aluminum, also in units of g/cm 2 , and E\ represents discrete energy values in 
units of MeV. These discrete source term values are used in the following way. Consider first 
the solution of equation (11) by the multigroup method for an all-aluminum shield with no 
target material; i.e., target thickness yj, — 0. The HZETRN program was run to simulate the 
solar particle event of February 23, 1956, and the source term .g{Ei,Xj,yk) associated with an 
aluminum- water shield was generated for these conditions. With this source term, equation (11) 
was solved by the multigroup method. 

For a single shield material, j3 = 1, equation (11) becomes 




fE fa\ 

<Kx,E) = J <7 S1 (£,£') <j>(x,E r ) dE' + g(E, x) 


( 36 ) 


where an integration of equation (36) from E; to Ej + \ produces 



"^dE + 

ox 


/■*(+ 1 

/ <r{E) <p(x,E) 

JEi 

rEj+i rE/u\ 

JEi JE 


dE 

{E,E7) <f>{x,E') dE' dE + 



dE 


( 37 ) 


We define the quantities 



( 38 ) 
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and interchange the order of integration of the double integral terms in equation (37). Then 
apply the mean value theorem to obtain the result 

d<t> f E '+' f E ' . , , 

+ v<t>i = / / (E,Ef) dE <j>(x,E ! ) dE’ 

dr JEi J E=E i 

f E i+2 f E i + 1 , , , 

+ / / a^{E,E')dE <p{x,E')dE l ' + b t (39) 

JE i+] J E =a ] E 1 


over the energy group Ej < E 1 < E; + j . For the energy spacing E i+ j = E-J a , the first double 
integral in equation (39) represents integration over the lower triangle illust rated in figure 2(c). 
The second double integral in equation (39) represents integration over the upper triangle 
illustrated in figure 2(c). Define 



cr S{ (E,E l )dE 
<r S] (E,E')dE 

E' 




(40) 


and then employ another application of a mean value theorem for integrals to write equation (39) 
in the form 

^ + TT*i = Sl [Ei + 9i(E i+ 1 - Ei)}<t>; + 9i [E i+ 1 + 0. 2 (E i+ 2 - E /+ i)]* i+ , + hi (41) 
ax 

This produces the coefficients associated with the energy group E t to £7+ i, which are given by 


flii = <71 “ 

+ 1“ 02 J 


In tins way, the diagonal and off-diagonal elements of the coefficient matrix in equation (26) are 
calculated. 

For a compound target material, comprised of material 1 and material 2, there are t wo values 
of a. A value a\ is selected for material 1 and a value is selected for material 2 of the 
compound material. In this case, equation (36) takes on the form 




fXjfa i 

o(x,E) = J (t^E^E 1 ) <j>{x,E') dE' 

.fE/a -2 

+ / <rs. 2 (E,E’) Hx,E')dE'+ g(E,x) 

J E 


(43) 


where cr A] and <r^ are scattering terms associated with the respective materials. These terms 
are calculated in the IIZETRN code. Two cases are considered. The first- case requires that 
the E/a -2 fine be above the E/a i line. (See fig. 2(d).) The second case is where «2 = 0 (the 
hydrogen case) and the limits of integration for the second integral goes to infinity. Each case 
is considered separately. 

For the first case, assume that a\ > «2 > 0 and select the exact energy spacing dictated 
by the E/a -2 line. Then proceed as for the single shield material. Integrate equation (43) 
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from Ei to E\- j_i and interchange the order of integration on the double integral terms. Define 

E- 

b t — y(E,x)dE and obtain the equation 

d<&; 

-j— -f <7$i = 1 11 + 1 12 + 7 21 + h'2 + *<’ 

where now 7 the 7-21 and I 22 integrals have, because of the exact spacings, the forms 

r£i+ 1 fE' 


(44) 


hi = f ’ +] / <t» 2 {E,E') dEcj>{x,E) dE 

JEj JE=E i 

r E i + 2 / ,£, J+1 , , . 

1-22= / <r, 2 (^) dE dE 

J Ei-L\ J E—c\‘)E } 


(45) 


Defining the terms 




f E' 

0) = / Vsj{E,E') dE 

J E—Ej 

O' = 1,2) 

0)=f‘ +l c t .(E,Ef)dE 

J E =o 7 E' 

0 = 1,2) 


and using the mean value theorem for integrals gives from equations (45) 

h\ — h\{'2)[Ei + @l(Ei+l ~ Ei)]$i 

and 

7 22 = *2(2)[ £, /+l + ^2(^+2 ~ £;+l)] $ (+l 

where 0 \ and 62 define intermediate energy values associated with the mean value theorem. 
The integrals and I\ 2 are associated with integration limits (E y E/a\) and energy intervals 
dictated by the selection of 02 for determining the exact energy spacings. These integrals are 
associated with the trapezoidal area 1 (Al) and triangular area 2 (A2) illustrated in figure 2(d). 
These areas are a fraction of the triangle areas associated with the line E 1 = E/a*}. These 
fractions are given by 


/l = 
/2 = 


0 + 1 - E ,) 2 - $(E i+l - Ej/cnHEj+i - o iC+l) ) 

— Ei+l-Ej)^ 

(Ej+i/ai - E; + \)(Ej + i - o iE; + i) 

(E;+l ~ Ej)(E i+ 2 ~ E j+l ) 


(46) 


and we write 


hi =h h i(i)$i 

7l2 = /2*2(l)$i+l 

The coefficients for the system of differential equations in equation (27) are then given by 

°11 =*1(2) + / l*l(l) _ & 

012 =*2(2) + /2*2(1) 


(47) 


( 48 ) 
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For the second case, of hydrogen, 02 equals 0; therefore one of the limits of integration be- 
comes infinite. Let aq determine the energy spacing in this case. Again integrate equations (45) 
over the energy interval (£’;,£',•+ 1 ), which is determined by the E' = E/a 1 line. With the def- 
initions given by equations (38), integrate equation (43) over the interval (Ej,Ej + \) and then 
interchange the order of integration in the resulting double integrals to obtain 


d<t>; 

— 


If + 1% + hi 


where 


I 


* 

1 


and 


n 


f E i + 1 fE rEj+2 rE i+ 1 

/ / a s JE,E')dE<j>(X'E , )dE'+ 

JEi J E=E; J E i+] J E=n\t 

' as l (E,E l ) dE <p(x,Ef) dE 1 

fE)+ 1 fE 1 [ E i+j+\ f E i 

/ / v n {E,E ) dE <j>(x,E ) dE' + T 

J Ej JEi fr\ ■> C'+j JEi 

+1 <t s . 2 (E,E) dE <j>(x ,E/)d Ef 


and for all N* greater than some integer N > 0, it is known that </>(x,E) will be zero. Define 


h(E?) = f J cr S] ( /•:,/•') dE 

JEj 

fE t+l 

h A (Ef) = / <t h {E,E?) dE 

J w] E 1 

fE' 

hrAEf) = J e <7s. 2 {E,E') dE 

h{j)= <r n (E,E)dE 

J Ei + j 


(E, < E' < £, +1 ) 
(E i+ i < E' < E i+ 2 ) 
(Ei<E’<E 1+ 1 ) 

(Ei+j < Ej < Ej+j+i) 


and then write the coefficients associated with the system of differentia] equations as 

a;,/ = h 3 + /15 - a 
a i,i + 1 = ^ l 4 + * 6 ( 1 ) 

<**,*+ 2= h (j(2) 
a i,i+ 3= ^0(3) 


a ij+n 


h 6(n) 


In this way a system of equations is generated that has the triangular form given by the 
system of equations in equation (27). 

Again use t he source term <i(Ej, x j , yj. ) obtained from the HZETR.N simulation of the 
solar particle event of February 23, 1956, associated with an aluminum-water shield-target 
configuration. Note that now the multigroup system of equations (eq. (27)) associated with 
equation (39) must be solved for the multiple atom target material of water. Consider the 
cases of discrete shield thickness x% 2 : 3 , ... and apply the multigroup method to the solution of 
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equation (11) applied to all target material y > 0. For each value of x/ considered, the initial 
conditions are obtained from the previous solutions generated where y — 0. This represents the 
application of the multigroup method to two different regions: region 1 of all shield material 
and region 2 of all target material. Then continue to apply the multigroup method to region 2 
for each discrete value of shield thickness, where the initial conditions on the start of the second 
region represents exit conditions from the shield region 1. This provides for continuity of the 
solutions for the fluence between the two regions. 

Results and Discussion 

The present formalism was used to evaluate the neutron fluence for various aluminum shield 
and water target combinations. Figure 7 illustrates the low-energy neutron fluence due to the 
scattering of evaporation neutrons in an aluminum shield for various thicknesses with y^ = 0 
(i.e., no target material). Figure 8 illustrates the total neutron fluence for various aluminum 
shield thicknesses. This fluence consists of the IIZETRN-generated neutron fluence plus the 
multigroup-generated low-energy neutron fluence. Figures 9, 10, and 11 are graphsof the neutron 
fluence in depths of 1, 10, and 100 g/cm 2 of aluminum generated from the HZETRN code both 
with and without the addition of the multigroup evaporation neutrons. 

Typical results for no shield before the water target are illustrated in figures 12, 13, and 14 
where a comparison of the multigroup method with the previous HZETRN results for thicknesses 
of 1, 10, and 30 g/cm 2 can be made. Note that in the calculations of the mult igroup method, the 
source terms g(E,x) y the scattering term t r^E.E / ), and cross section cr(E) of equation (11) are 
all given as numerical output from the HZETRN code for the solar particle event of February 23, 
1956. Also note that these calculations were compared with the LAITET Monte Carlo results 
from reference 24. Figures 12, 13, and 14 illustrate this comparison for neutron fluences versus 
energy at water depths of 1, 10, and 30 g/cm 2 , respectively. Figure 15 is a graph of neutron 
fluence versus depth in a shield-target configuration of 100 g/cm 2 of aluminum followed by 
100 g/cm 2 of water. Observe the increase in the low-energy neutron fluence at the aluminum- 
water boundary. This increase is caused by high-energy neutrons colliding with hydrogen atoms, 
which results in large energy losses. In these types of collisions, the neutrons of modest energies 
give up one half of their energy on the average; thus, the lower energy neutron fluence is increased. 

In figures 12, 13, and 14, note the distinct improvement of the fluence by using the multigroup 
evaporation neutrons over that of the previous HZETRN results. These improved results are 
still a little lower than the results predicted by the Monte Carlo simulation. These figures 
show that the multigroup method is more accurate at the higher target depths compared with 
results at the lower depths. This is due to the straight -ahead approximation assumptions 
used in the one- dimensional Boltzmann equation, where all secondaries produced by nuclear 
collisions are assumed to move in the same direction as the primary nucleon which caused the 
collision. This assumption is true for secondaries which are high-energy particles. This straight- 
ahead approximation is not true for low-energy neutrons produced by evaporation because these 
neutrons are generally isotropically distributed. These neutrons make up the source terms in the 
multigroup method. The straight-ahead assumption causes errors at the smaller target depths 
because it fails to account for all the low-energy neutrons transported back from larger depths of 
the material. In an attempt to improve the performance of the multigroup method for simulating 
low-energy neutrons, the assumption was made that only one half the source terms moved in 
the forward direction while the other half moved in the backward direction. The solution of the 
multigroup system of equations (eq. (27)) was then modified. Using one half the source terms 
g(Ej,Xj } yfc) y system of equations (eq. (27)) was marched first through the shield material and 
then through the target material. By using the end boundary condition generated, the equations 
were then marched backwards through the target and then the shield material. The fluences 
from the forward and backward marching were then added to obtain a total fluence. This process 
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is referred to in the figures as the two-dimensional multigroup method. Figures 16, 17, and 18 
illustrate the results of the two- directional multigroup method applied to the case of no shield 
and a target of water only for nominal depths for an exposure to the solar particle event of 
Februaiy 23, 1956. Figure 19 illustrates the fluence in a depth of 10 g/cm 2 of water when the 
two-dimensional method was applied to a 100 g/cm 2 aluminum shield followed by a 100 g/cm 2 
target of water when exposed to the solar particle event, of February 23 , 1956. Observe that the 
two-directional multigroup method greatly improves the low-energy fluence predictions at the 
smaller depths. 

Research is continuing to dost' the remaining gap between transport code predictions and 
Monte Carlo results. Possible errors from various sources are being investigated. The nuclear 
cross sections used are believed to be one source of error because only elastic cross sections 
were used in the multigroup simulation. The elastic cross sections are much larger than the 
nonelastic cross sections at low energies. Nonelastic cascading does occur and it is believed 
that the mult igroup method would be improved by incorporating both types of cross sections. 
Other sources of errors reside in the HZETRN program itself. The nuclear cross sections used 
by HZETRN arc interpolated from a large database that was developed experimentally many 
years ago; this database needs to be updated. The HZETRN code is a one-dimensional transport 
code using the straight-ahead approximation. The improvement of the multigroup method in 
going from the straight-ahead approximation to the two-directional multigroup approximation 
suggests that similar type changes be incorporated into the HZETRN code in order to reflect 
the nonisotropic character of the events. 

Concluding Remarks 

These preliminary studies have shown that the multigroup method developed for the study of 
low-energy neutron transport has made significant improvements in and is compatible with the 
current HZETRN code developed at Langley Research Center. It has proven to be a fast and 
efficient algorithm for the inclusion of low-energy neutrons into the HZETRN code. The addition 
of nonelastic processes in the low-energy neutron transport is expected to further improve the 
result. 
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Figure 7. Low-energy neutron fluences due to scattering of evaporation neutrons in aluminum 
shield exposed to solar particle event of February 23, 1956. 



Energy, MeV 


Figure 8. Total neutron fluences in aluminum shield exposed to solar particle event of 
February 23, 1956. 





Figure 9. Energy spectra of neutron fluence at depth of 1 g/ cm 2 in aluminum exposed to solar 
particle event, of February 23, 1956. 



Figure 10. Energy spectra of neutron fluence at depth of 10 g/cm 2 in aluminum exposed to solar 
particle event, of February 23, 1956. 
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Energy, MeV 


Figure 11. Energy spectra of neutron fluence at depth of 100 g/cm 2 in aluminum exposed to 
solar particle event of February 23, 1956. 



Figure 12. Energy spectra of neutron fluence at depth of 1 g/cm 2 in water exposed to solar 
particle event of February 23, 1956. 
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Figure 13. Energy spectra, of neutron fluence at depth of 10 g/cm 2 in water exposed to solar 
particle event of February 23, 1956. 
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Figure 15. Neutron fluence as function of depth in 100 g/cm 2 aluminum shield followed by 
100 g/cm 2 water target. 



Figure 16. Energy spectra of neutron fluence at depth of 1 g/cm 2 in water exposed to solar 
particle event of February 23, 1956, and calculated with two-directional multigroup method. 




Figure 17. Energy spectra of neutron flue nee at depth of 10 g/cm 2 in water exposed to solar 
particle event of February 23, 1956, and calculated with two-directional multigroup method. 



Figure 18. Energy spectra of neutron fluence at depth of 30 g/cm 2 in water exposed to solar 
particle event of February 23, 1956, and calculated with two-directional multigroup method. 
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Figure 19. Energy spectra of neutron fluence at depth of 10 g/cm^ in target of 100 g/cm 2 
aluminum shield with target of 100 g/cm 2 of water behind it when exposed to solar particle 
event of February 23, 1956, and calculated with two-directional multigroup method. 
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